The invention relates to seismic exploration for oil, gas and other underground resources and to improving seismic mapping of subsurface formations by the use of a deconvolution technique which gives good results even when the observed seismic signals are not minimum phase or highly non-Gaussian. In a broader sense, one aspect of the invention relates to improving a signal resulting from the physical interaction (convolution) of an unknown wavelet with an unknown stationary random sequence such that the smoothing effect of the wavelet can be removed, or at least substantially removed, and the random sequence can be revealed. In that sense, the wavelet resulting from the explosion which can be used in seismic exploration is an example of an unknown wavelet, the positions of the subsurface seismic reflectors are an example of the members of a stationary random sequence, and the signal to be improved is the seismic trace based on surface measurements of the seismic signal resulting from the interaction (convolution) of the wavelet caused by the explosion and the subsurface reflectors (at the interfaces of subsurface layers having different acoustic properties). Other aspects of the invention are addressed in the detailed description below.
Subsurface reflection seismology involves mapping of subsurface formations on the basis of the arrival times of events reflected from subsurface layers. Seismic energy generated at the surface penetrates the earth's layered media and some of it is reflected back at the interfaces between layers. The reflected energy arriving at the surface is measured by receivers (detectors at the surface. On land, seismic energy typically is generated by chemical explosions, vibrating machines or thumping devices. Preplaced receivers, also called detectors or geophones, arranged in an array or "spread," detect the seismic energy which comes back. At sea, a source such as an array of air guns is actuated every few seconds as the ship moves over a predetermined course. The returning seismic energy is picked up by detectors embedded in a cable (a streamer), trailing the ship.
The complexity of the subsurface structure and the noise and other imperfections of the observed signals make it necessary to process these signals through a number of operations, each involving serious approximations and shortcuts, before a useful map of the subsurface formations can be deduced. One major operation in such signal processing is velocity analysis, made necessary by the fact that the seismic signal observed at the surface shows the reflection events as a function of time, and must be converted into a function of depth in order to derive a more useful map of the subsurface structure. However, seismic waves have a velocity which is very much dependent on the propagation medium, and velocity changes significantly as the waves travel into the earth. Generally, the velocity increases with depth, although occasionally there may be layers in which a decrease in velocity occurs. For a given surface point, the velocity plotted as a function of depth is typically called the velocity function. Thus, there are two important variables: time of reflected events and velocity. From knowledge of these variables, the depth to the reflecting horizons can be determined. Because there are important lateral changes in velocity, and the velocity function varies from one location to another, a given velocity function cannot be assumed to be valid for the entire prospect but must be continuously corrected from place to place over the area of exploration. When a deep enough borehole is available in the right place and can be logged, the velocity function can be estimated from the outputs of seismic detectors placed in the hole at various depths, or by vertical seismic profiling ("VSP"). In most cases, however, it is necessary to estimate the velocity function by measurements confined to the surface, and by subjecting these measurements to sophisticated signal processing. Velocity analysis results are also used in the so-called dynamic corrections outlined immediately below.
Another major operation in processing the observed signals involves corrections applied to the receiver outputs. The so-called static corrections include a source correction and a detector correction, the combined effect of which is to put the seismic signal source and the seismic signal detectors in a fictitious horizontal plane. In practice, this is done by appropriately time shifting the respective receiver outputs. The dynamic corrections convert each receiver output to the output that would have been received if the source and receiver were at the same lateral point, e.g., the point midway between the actual source and receiver locations. In this conversion the traces are considered as made up of primary reflections, i.e., reflections resulting from a raypath from the source down to the reflecting horizon and then directly up to the receiver. If the layers are horizontal, then all the reflection points (or depth points) are directly beneath the midpoint between source and receiver. If the layers are dipping, then the depth points are offset from this midpoint. The dynamic corrections therefore depend on both the velocity function and the dip of the reflecting beds. The component of correction due to the separation of source and receiver is called the normal moveout correction, and the component due to dip is called the dip correction. Thus, four important corrections to the recorded traces (detector outputs) are the source corrections, receiver corrections, normal moveout corrections, and dip corrections.
Each observed signal (recorded trace) can be considered as a time series made up of reflected events together with various interfering waves and noise. The desired reflected events are the primary reflections, that is waves that travel directly down to a given reflector and then directly back up to the surface where they are recorded. An important type of undesired interfering wave is the multiple reflection, which has a raypath that goes down to a given reflector, then up to another reflector, then down to a reflector, then up again, possibly with additional up and down paths, until the multiply reflected event reaches the surface and is recorded. In any layered system there can be infinitely many possible multiple reflections. The presence of such multiple reflections makes the identification of primary reflections difficult. Accordingly, it is necessary to attenuate the multiples as much as possible.
To facilitate the extraction of the useful signal, typically multifold coverage is recorded by laying out a source and a spread of detectors, activating the source and recording the detector outputs, then moving the entire configuration laterally, along the detector line, and repeating the process. When the configuration is moved in small enough increments between shots, each depth point can be covered several times, thereby introducing considerable redundancy which allows subsequent enhancement of the primary events and suppression of multiple events. If all the traces in the prospect are segregated or gathered into sets called gathers such that all traces within a given gather have a common midpoint between source and receiver, the appropriate moveout correction should convert each trace in a given gather into about the same equivalent trace, namely that primary reflection trace which would have been received if source and receiver were directly at the common midpoint in question. Stated differently, the primary reflections of all the traces after appropriate moveout should tend to be in phase while the multiples should tend to differ as between traces and therefore should typically be out of phase. If this approach is used for the normal moveout correction first, then it can also be applied to making the appropriate source, receiver and dip corrections by different gatherings of the traces. Alternately, a simultaneous corrections approach can be used, either directly or through using iterative techniques. If the corrected traces in a common midpoint gather are added up, the result should be severe attenuation of the multiple events, because they should be out of phase with each other, and amplification of the primary events, which should be in phase with each other if the gathers were properly selected and corrected. Thus one output trace for each midpoint can be produced which is called the stacked trace for that midpoint.
A special kind of multiple reflection is the so-called reverberation. One conceptualization of reverberation is energy which gets trapped in the near surface layers and keeps being reflected up and down. This energy can become attached to the primary reflections as they travel through the near surface layers. As a result, instead of being sharp and clear and with good time resolution, the observed reflections are followed by long reverberation trains. These trains overlap with the trains of succeeding reflections, and thus the whole seismic trace (observed signal) can be given a ringing or sinusoidal character which makes it difficult or impossible to pick out the onset of the primary reflections. Cancelling the energy of each reverberation train, but leaving intact the primary reflection, can increase the resolution of all the reflected events. This process of cancelling the reverberation trains typically is called deconvolution, and can be thought of as the inverse of the geophysical process of "convolving" (i.e., rolling together or interacting) the source wavelet with the subsurface formations. However, as will become clearer below, this is not the mathematical process of deconvolution.
One known type of deconvolution is based on considering a wavelet of the observed signal as made up of a primary reflection with its attached reverberation. With this assumption, the wavelet should be minimum delay, in the sense that most of its energy is concentrated in its initial part. Moreover, all such wavelets should be similar in overall shape. Because the primaries result from geologic beds believed laid down with irregular thicknesses, the arrival times of the primary events should be effectively random. Hence, the auto-correlation function of the trace should be about the same as the auto-correlation function of the wavelet, and so from this auto-correlation function one can estimate the required inverse (deconvolution) operator. The application of this operator to the trace should yield the deconvolved trace, namely a trace where the reverberation components of the wavelets have been removed, thereby increasing the resolution of the primary reflections. The process of deconvolution can be extended to remove certain long period multiple reflections as well. In practice, deconvolution can be applied to the seismic traces either before or after stacking, depending upon cost and other considerations.
After deconvolution, the processed signals typically are subjected to a technique called migration or depropagation, which can be conceptualized as a process of propagating the waves observed on the surface backward in time into the earth to the underground structure, to reveal the actual spatial position of the subsurface reflection points at depth. This process can also be described as the transformation of signals observed at the surface to signals which would have been observed at depth.
As discussed in an article by the inventor herein (Wiggins, R. A., Minimum Entropy Deconvolution, Geoexploration, 16 (1978) 21-35), if the source wavelet is known at least approximately, then wavelet shaping by estimating a least squares or Wiener inverse of the source wavelet can lead directly to a good estimate of the desired signal. If the source wavelet is minimum phase or minimum delay, and the desired signal is approximately white, then predictive deconvolution can be used. If the amplitude spectrum of the source wavelet is smooth with respect to that of the desired signal, then homomorphic deconvolution can be utilized. When the desired signal consists of a few large spikes and the ensemble of sample signals exists for which the time separations between the spikes differ while the source wavelet shape remains constant, then minimum entropy deconvolution can be used in the manner described in detail in the article. Further background information can be found in the references cited at pages 34 and 35 in the Wiggins article. Yet further background information can be found in Robinson, E. A., et al., Deconvolution of Geophysical Times Series in the Exploration for Oil and Natural Gas, Elsevier Scientific Publishing Co., 1969 and in Robinson, E. A., Time Series Analysis and Application, Goose Pond Press, 1981. The Wiggins article and said Robinson and Robinson et al. books are hereby incorporated by reference in this specification.
The known "deconvolution" processes applied to seismic traces are believed to have a number of shortcomings. For example, they are not believed to work well with seismic traces which are not minimum phase (delay) or sufficiently non-Gaussian, and different signal processing techniques must be applied depending on which one is the case. Moreover, it is believed that the known processes suffer from complexity and insufficient accuracy.
Accordingly, an important object of this invention is to provide an improved technique of deconvolving the result of the convolution of a wavelet and a non-Gaussian random process, such as the result of the interaction of a seismic source wavelet with the subsurface reflectors, in a manner which works well for all or most typical seismic traces, uses the same approach regardless of the type of trace, processes a trace only as long as needed to deconvolve it sufficiently, and operates rapidly and efficiently. In a broader sense, an important object is to find a practical solution to the so-called blind deconvolution problem, encountered when an observed or measured signal can be considered for practical purposes to be the result of the convolution (interaction) of an unknown wavelet with a stationary random sequence. The goal is to remove the effects of the wavelet to an extent which would allow the underlying random sequence to be revealed and examined free, or substantially free, from the smoothing effects of the wavelet. The invention is particularly useful in solving in a practical manner a long standing problem in seismic exploration, but is also useful in other fields in which the observed signal is not the one of interest, but is the convolution of a wavelet with the sequence of interest, such as, without limitation, in certain sonar and radar activities and in other fields making use of acoustic signals and other signals which get contaminated by an interaction with wavelets which can be treated for practical purposes as convolution.
In conceptual terms, and as applied to seismic exploration, the improved deconvolution technique involves finding what can be characterized as an inverse filter which, when applied to the observed signal, changes it to a good estimate of the reflectivities and locations (in time) of the depth points producing primary reflections giving rise to the observed signal. In practicing one nonlimiting example of the invention, this is done by successively modifying the observed signal by two-term operators. The first operator is chosen on the basis of the observed trace and beliefs as to the entropy of the expected reflectivity series, and is applied to the observed signal to produce a modified observed trace. The second operator is chosen on the basis of the modified observed trace and similar entropy considerations, and is applied to that modified trace to produce a new modified trace. The third operator is chosen on the basis of the last modified observed signal and similar entropy considerations and is applied to the last modified observed signal to produce yet another modified observed signal, and so on until no significant further improvement is apparent. The latest modified observed trace is then a good estimate of the reflectivities and locations (in time) of the subsurface reflectors of interest. Each two-term operator can be of the form [1,0,0, . . . ,a] or [a,0,0, . . . ,1], where the separation between its 1 and a terms is determined from the lag of the maximum deviation of the auto-correlation of the most recently modified (if at all) observed signal, and the amplitude of the term a is found from other characteristics of the most recently modified observed signal (e.g., its auto-correlation). Which of the possible operators is used in a given run through the process is determined on the basis of which gives less entropy of the modified observed signal, and this can be estimated, for example, by using a Varimax measure.
The major steps of a particular nonlimiting example of practicing the invention, as applied to seismic exploration, are as follows:
1. Generate an observed signal. If a stacked trace is used for the observed signal, in the first pass through step 1 the observed signal is the stacked trace itself, and in each subsequent pass through this step 1 the observed signal is the result of most recent pass through step 5 below. The observed signal can be the digitized signal y(i), where i=1,2, . . . ,I, and the index i designates digital samples taken over a period long enough to include the deepest reflector of interest (typically several seconds) and the sampling frequency is sufficiently high to satisfy the Nyquist criteria for the highest observed signal frequency of interest. While for simplicity only the example of using the stacked trace as the observed signal will be discussed in detail, this need not be the case. The observed signal can be, as nonlimiting examples, the output of one receiver, some combination of receiver outputs other than the stacked trace, and it can be an observed signal in surface reflection seismics or the observed signal in other types of seismic exploration, such as vertical seismic profiling (VSP), or it can be some other measured signal, such as a sonar or a radar signal;
2. Find the auto-correlation function R(y,y;i) of the observed signal y(i) for i=0,1,2, . . . ,I;
3. From the auto-correlation R found in step 2, find the lag d of a selected side lobe of the auto-correlation (e.g., the distance, expressed in the number of steps i, between the central lobe of the auto-correlation, at R(y,y;0) and a selected side lobe, at R(y,y;i). This lag d determines the number of zeroes between the two terms of the operator (if d=i, then there are (d-1) zeroes);
4. From the same auto-correlation R(yy) found in step 2, find the amplitude of the term a. One example of a relationship useful for finding the amplitude of the term a is a=[-R(yy;d)]/[R(yy;0)]. R(yy;d) is the amplitude of the auto-correlation function of the current observed trace at the selected side lobe and R(yy;0) is the auto-correlation amplitude at the central lobe. Other relationships between auto-correlation attributes can be used instead, so long as the resulting term a, when used to operate on the current observed signal as the first or last term of a two-term operator, minimizes (or at least significantly reduces) the entropy of the result as determined, for example, by checking if a Varimax measure V is maximized (e.g., by checking if the new V is greater than the previous one plus a threshold, and continuing the process only if the answer is affirmative). Having determined the separation d (and, hence, the number of zeroes between 1 and a), and the amplitude of a, one of two possible two-term operators [1,0,0, . . . ,a] and [a,0,0, . . . ,1] is selected on the basis of which is more effective in reducing the entropy of the observed signal resulting from its application to the current version of the observed signal. This can be done, for example, on the basis of a Varimax measure V such as V=y(i).sup.4 /[y(i).sup.2 ].sup. 2, where y(i) is the modified observed signal resulting from the application of the two-term operator being tested to the current version of the observed signal. The two-term operator giving the larger V is selected, as it is the one determined to be more effective in reducing the entropy of the observed signal; and
5. Apply the two-term operator chosen in steps 3 and 4 to the observed signal to generate a modified observed signal. In the first pass through this step 5 the two-term operator is applied to the actual observed signal, e.g. a stacked trace; in each subsequent pass through this step 5 the newly selected two-term operator is applied to the previous result of this step 5.